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Abstract 

Moving  tracked  vehicles  excite  large-amplitude  seismic  surface  waves  that  can  be  used  to  track  and 
identify  them  at  ranges  over  1  km.  Furthermore,  these  surface  waves  generally  possess  robust  spatial  coherence, 
show  a  smooth  amplitude  decay  as  a  function  vehicle  range,  and  are  minimally  affected  by  severe  meteorological 
conditions.  Because  of  these  properties,  seismic  signals  should  be  used  to  augment  acoustic  sensing  in  battlefield 
systems.  However,  large  changes  in  vehicle  signature  characteristics  can  be  produced  by  geological  variations.  The 
heightened  interest  in  using  seismic  signals  for  battlefield  applications  has  created  a  need  to  understand  the  complex 
effects  produced  by  the  ground  on  propagating  seismic  surface  waves.  High  fidelity  forward  modeling  can  be  used 
to  both  explain  these  effects  and  to  provide  raw  data  for  system  development.  Using  synthetic  data  in  this  manner 
can  reduce  system  development  time  and  overall  costs,  while  simultaneously  improving  system  performance. 

Geologic  inhomogeneity  and  material  properties  affect  signal  characteristics  at  target  ranges  as  short  as 
100  m;  signals  from  more  distant  target  ranges  are  affected  to  an  even  greater  extent.  Horizontal  inhomegenities 
resulting  from  subsurface  variation  and  topography,  in  particular,  affect  seismic  surface  wave  characteristics.  The 
need  to  accommodate  strong  near-surface  inhomogeneity  and  seismic  body-wave  conversions  at  geologic 
boundaries  compels  the  use  of  discrete  numerical  propagation  approaches.  A  finite  difference  time  domain  (FD-TD) 
approach  was  selected  as  the  propagation  method  for  simulating  seismic  signatures  because  of  the  wide  availability 
of  sophisticated  FD-TD  codes  and  their  complete  consideration  of  all  3D  body  and  surface  wave  energy  transfer 
processes. 

We  discuss  the  mathematical  basis  of  our  FD-TD  code  and  present  simulated  propagation  results  from  a 
large  parallel  3D  FD-TD  elastic  model.  For  portability  reasons,  the  code  is  written  in  standard  Fortran77  and  is 
parallelized  using  the  Message  Passing  Interface  (MPI)  subroutine  library.  The  code  has  been  successfully  run  on 
Sun  workstation  clusters  and  on  massively  parallel  Cray  and  IBM  platforms.  Simulations  are  discussed  for  a  plane 
layered  geology  and  for  a  geology  with  dominant  3D  features.  Both  results  use  an  explosive  pressure  impulse  placed 
just  under  the  earth’s  surface.  The  3D  geologic  model  contains  an  isolated  hard  rock  topographic  feature  protruding 
through  a  layered  near-surface  soil.  In  both  cases  complex  seismic  waves  are  observed.  In  the  case  of  the  protruding 
topographic  feature,  the  results  show  strong  surface  wave  reflection  and  refraction  around  it.  These  phenomena 
dramatically  affect  the  character  of  seismic  signals  by  altering  the  waveform,  by  reducing  the  signal  amplitude,  and 
by  reducing  the  spatial  coherence  of  the  wavefields.  Signal  effects  of  this  magnitude  would  severely  impact  system 
performance.  Conversely,  foreknowledge  of  these  effects  can  be  used  advantageously  to  optimally  place 
autonomous  battlefield  monitoring  systems. 


1.0  Introduction 


Use  of  synthetic  (modeled)  data  can  substantially  reduce  the  time  and  costs  of  developing  a  combat 
system,  while  simultaneously  improving  its  reliability.  This  is  achieved  by  considering  how  the  system  will  perform 
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under  diverse  deployment  scenarios  and  by  discovering  physical  processes  that  exhibit  robust  environmentally 
invariant  signal  characteristics  that  are  specific  to  a  target  class. 

The  central  long  term  objective  of  our  project  is  to  produce  high  fidelity  seismic  signatures  generated  by 
heavy  ground  vehicles,  such  as  armor,  motorized  artillery,  or  TEL  transports,  operating  in  complex  geologic  settings 
such  as  that  shown  in  Figure  1.  The  approach  must  be  flexible  enough  to  accommodate  a  wide  variety  of  vehicles, 
background  noise,  and  geologic  settings.  Simulations  should  also  consider  target  ranges  of  approximately  1  km  or 
more.  Greenfield  and  Moran  (1999),  Moran  et  al.  (1997,  1998),  and  Prado  (1998)  have  shown  that  the  seismic 
signatures  of  tracked  ground  vehicles  are  dominated  by  surface  waves  that  often  contains  multiple  modal  arrivals. 
Furthermore,  the  target  vehicles  of  interest  are  spatially  large  in  extent  (relative  to  a  wavelength)  and  in  general  they 
apply  both  vertical  and  horizontal  tractional  and  pressure  forces  to  the  earth  or  road  surface. 

Generating  a  synthetic  seismic  signature  requires  specification  of:  1)  vehicle  forcing  functions  on  the  earth, 
2)  definition  of  the  local  geology,  and  3)  propagation  of  seismic  energy  away  from  the  vicinity  of  the  target  to  a 
monitoring  sensor.  Each  of  these  major  simulation  components  presents  substantial  technical  challenges.  In  this 
paper,  we  report  on  a  FD-TD  method  of  seismic  energy  propagation. 

Elastic  materials  readily  support  multiple  seismic  wave  types,  including  compressional  (P)  and  shear  (S) 
body  waves,  as  well  as  a  variety  of  modal  surface  wave  phases.  The  later  typically  are  strongly  dispersed.  Geologic 
material  properties  dramatically  affect  signal  characteristics  at  target  ranges  as  short  as  100  m;  signals  from  more 
distant  target  ranges  are  affected  to  an  even  greater  extent.  Horizontal  and  vertical  geologic  inhomegenities  resulting 
from  depositional  layering  and  topography  affect  seismic  surface  wave  characteristics  by  creating  lossy  waveguides, 
by  scattering  signals,  and  through  conversions  of  incident  wave  phases.  The  need  to  accommodate  3D  geological 
inhomogeneity  and  complex  wave  conversions  compels  the  use  of  a  discrete  numerical  propagation  approach.  A 
variety  of  numerical  methods  have  been  reported  in  the  geophysical  literature.  These  include  wavenumber 
integration  (only  applicable  for  plane  layered  geologies),  finite  element  (FE),  boundary  element  (BEM),  and  FD-TD 
techniques.  The  wide  availability  of  expertise  and  the  comparative  sophistication  and  maturity  of  FD-TD  codes 
make  this  approach  the  most  attractive  numerical  solution  method  for  our  problem. 


1  Km 


Figure  1.  A  moving  armored  vehicle  generates  large-amplitude 
seismic  signals.  Complex  geology  has  a  strong  effect  on  signal 
characteristics.  This  study  applies  3D  FD-TD  propagation  to  consider 
the  dynamic  source  characteristics,  complex  geology,  and  spatially 
distributed  noise  fields. 


2.0  Solution  Method 

We  provide  a  brief  summary  of  the  theoretical  and  numerical  methods  employed  to  solve  the  propagation 
problem.  Greater  detail  can  be  found  in  Hestholm  and  Ruud  (1998).  We  implement  a  3D  FD-TD  elastic  wave 
equation  model  using  a  staggered  grid,  central  difference  scheme.  Propagating  fields  are  expressed  with  a  velocity- 
stress  parameterization.  The  differential  operator  expansions  are  8th  order  accurate  in  space  and  2nd  order  accurate  in 
time.  A  curvilinear  coordinate  system  is  used  to  incorporate  the  effects  of  surface  topography  on  seismic  energy 
transfer  processes.  On  the  surface  of  the  earth  (called  the  free  surface),  we  impose  a  zero  stress  boundary  condition. 
The  remaining  gird  boundaries  are  terminated  using  an  exponential  decay  condition. 


2.1  Equations  of  Motion 


To  consider  the  effects  of  topography,  we  incorporate  Hestholm  and  Ruud’s  (1998)  curvilinear  coordinate 
system  transformation.  The  approach  begins  by  specifying  the  surface  of  the  earth  as  an  arbitrary  single-valued 
elevation  function  z 0(x,y).  Consider  a  curved  system  with  spatial  coordinates  given  by  £  r,  and  77.  The 
transformation  is  defined  as 


£  =  X,  1. 

T=y,  2. 

V  =  ^—z0(x,y),  3. 

^  Max 

where  Tj  is  the  depth  in  the  curved  coordinate  system,  z  is  the  depth  below  the  surface  in  the  rectangular  system,  and 
Zmux  is  the  maximum  depth  of  interest.  The  transform  is  illustrated  in  Figure  2.  It  effectively  stretches  or  shortens  the 
vertical  coordinate  (77)  corresponding  to  the  variation  in  z0-  Using  to  represent  a  continuous  wave  field  in 

the  curved  coordinate  system,  we  apply  equations  1-3  with  the  chain  rule,  which  leads  to 
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Using  1  -  3  we  define, 
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Applying  equations  1  -  9  to  the  3D  inhomogeneous  elastic  wave  equation  leads  to  the  following  particle 
velocity  expressions 
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where  u.  v,  and  w,  are  the  vector  particle  velocities,  p  is  the  material  density,  are  body  forces,  and  Gy  are  the  stress 
fields.  For  isotropic  propagation,  the  3-by-3  stress  matrix  is  symmetrical  around  the  diagonal  and  has  only  6  unique 
components.  These  are 
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where  ju  is  the  shear  modulus  and  A,  Lame’s  constant.  Equations  10-18  completely  specify  all  seismic  wave 
phases  (compressional,  shear,  and  surface  waves)  for  an  isotropic,  purely  elastic  system.  The  equations  of  motion 
given  above  do  not  consider  energy  loss  through  friction,  or  visco-elastic  response.  For  propagation  in  soils,  these 
are  important  effects  that  will  be  included  in  later  investigations  through  further  generalization  of  equations  10-18. 


Geologic  Model 


Computational  Domain 


Figure  2.  Curvilinear  coordinate  system  transform.  A)  The  geological  model  and  the  propagating  wave  fields 
are  specified  in  a  curved  system  whose  upper  surface  conforms  with  local  topography.  The  grid  is  pinched  or 
dilated  in  the  vertical  dimension.  B)  Propagating  wave  fields  are  transformed  into  a  rectangular  system  and 
standard  FD-TD  expansions  are  applied. 


In  addition  to  the  traditional  time  stability  criteria  requirements  (Lines  et  al.,  1999)  of  FD-TD  methods,  the 
curvilinear  coordinate  basis  of  this  solution  approach  imposes  a  spatial  stability  requirement.  It  can  be  readily  seen 
in  equations  10-18  that  if  the  constants  A,  B ,  and  C  are  larger  than  1,  as  time  progresses,  the  solutions  will  show 
roughly  geometric  growth.  To  ensure  stability,  the  slopes  of  the  bounding  topographic  surface  should  be  small 
relative  to  the  thickness  of  the  model  space.  This  condition  is  met  when 


The  Lame  material  property  parameters  are  related  to  the  compressional  (P)  and  shear  (S)  body  wave 
propagation  speeds  (V)  by 


Density,  VP,  and  Vs  are  readily  observed  through  direct  measurement.  Thus,  these  parameters  are  used  as  input  to  the 
geologic  model. 

In  our  numerical  code,  the  differential  operators  in  equations  10  -  19  are  expressed  via  discrete  8th  order 
finite  difference  expansions  in  space  and  2nd  order  expansions  in  time.  For  the  purpose  of  portability,  the  code  is 
written  in  Fortran77.  FD-TD  methods  require  large  computation  resources.  To  make  the  numerical  execution 
practical,  the  algorithm  has  been  parallelized  using  a  domain  decomposition  strategy.  This  approach  apportions  the 
numerical  model  into  a  number  of  computational  sub-domains  that  are  associated  with  a  dedicated  CPU.  As  shown 
in  Figure  3,  a  single  processor  is  usually  assigned  to  each  sub-domain.  This  effectively  divides  the  size  of  the 
computation  by  the  number  of  processors.  Exchange  of  data  among  processors  is  accomplished  by  subroutine  calls 
using  the  Message  Passing  Interface  (MPI)  library.  MPI  is  in  the  public  domain  and  has  become  a  popular  standard 
for  parallel  numerical  problems.  One  of  its  most  appealing  features  is  its  platform  independence.  CRREL’s  parallel 
code  has  been  successfully  run  on  Sun  workstation  clusters,  and  on  massively  parallel  Cray  T3E  (512  CPU  nodes) 
and  IBM  SP2  (256  CPU  nodes)  platforms. 


A)  Full  Model  Domain 


B)  Decomposed  Domain 
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Figure  3.  A)  A  2D  Model  on  a  12-by-12  rectangular  grid.  B)  The  model  decomposed  into  9  sub-domains. 
Each  sub-domain  spans  a  4-by-4  grid.  Each  of  these  9  domains  is  assigned  to  a  single  processor  that  is 
situated  in  a  3-by-3  processor  mesh  topology.  MPI  provides  protocols  for  exchanging  data  among  processors. 
In  real  cases  typical  models  span  millions  of  grid  nodes  that  are  decomposed  and  evaluated  on  tens  of 
processors. 


2.2  Boundary  Conditions 

FD-TD  calculations  support  wave  propagation  within  a  bounded  region.  The  bottom  and  lateral  grid 
terminations  are  generally  not  physically  real  and  as  a  consequence  they  introduce  unwanted  reflections.  The  upper 
boundary  of  the  computational  grid  is  defined  as  the  earth  air  interface.  This  is  a  naturally  occurring  termination  and 
is  handled  by  the  imposition  of  free-surface  boundary  conditions. 

In  elastic  calculations,  the  surface  of  the  earth  represents  a  comparative  vacuum  that  requires  the 
imposition  of  special  boundary  conditions.  This  is  an  important  aspect  of  the  simulation  problem,  since  the 
imposition  of  these  special  conditions  is  largely  responsible  for  the  occurrence  and  characteristics  of  seismic  surface 
waves.  Stresses  can  not  be  supported  in  a  vacuum.  In  the  curved  coordinate  system,  this  requires 
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To  apply  these  stress  release  conditions,  the  normal  at  each  surface  grid  point  is  calculated  from  the  local 
topography  and  the  coordinate  system  transformation  is  applied.  Details  of  this  rotation  and  transformation  are  given 
by  Hestholm  and  Ruud  (1998).  The  result  is  a  system  of  linear  equations,  which  are  solved  by  direct  matrix 
inversion  using  2nd  order  vertical  and  horizontal  particle  velocity  derivatives  from  the  grid  nodes  immediately  below 
the  free-surface. 

Because  of  practical  limitations  on  computational  resources,  the  spatial  extent  of  the  bounded  region  is 
often  undesirably  small.  In  the  confined  space,  the  seismic  wavefronts  often  reflect  from  the  terminated  grid  with 
comparatively  large  amplitudes  that  do  not  correspond  to  reflections  observed  in  the  physical  world.  There  are  a 
variety  of  numerical  methods  that  can  be  applied  to  reduce  these  unwanted  boundary  reflections.  The  method 
implemented  here  follows  Cerjan  et  al.  (1985).  In  this  approach  an  exponential  reduction  is  applied  to  the  field 
values  inside  an  annular  shell  of  finite  width  ( NAx )  whose  outer  surface  is  the  edge  of  model  domain.  The  reduction 
is  given  by 

G(/i)  =  e“k»*g  25. 

Cerjan  et  al.  (1985)  recommends  a=0.015  and  A=20.  For  grid  positions  inside  the  annular  shell,  each  particle 
velocity  and  stress  field  component  is  scaled  by  G  at  every  time  step.  The  resulting  spatial  taper  is  very  slight  over 
N;  however,  the  reduction  occurs  at  each  time  step  for  both  the  inwardly  propagating  wave  and  its  returning 
reflection.  We  have  found  that  reflections  are  effectively  reduced.  However,  there  is  a  substantive  loss  in  usable 
computational  volume. 


3.0  Results 

A  model  run  progresses  for  a  user  specified  number  of  time  steps.  At  each  time  step  the  particle  velocity 
and  stress  fields  are  updated  at  every  grid  point  in  the  model  volume  and  are  thus  available  as  output.  Owing  to  the 
large  amount  of  data  involved,  it  is  frequently  desirable  to  save  only  a  subset  of  the  available  volumetric  space.  In 
this  paper  we  discuss  particle  velocity  results  from  vertical  or  horizontal  slices  taken  from  the  available  3D  model 
volume  at  selected  times.  These  are  termed  snapshots.  Corresponding  animated  results  are  discussed  in  the 
presentation.  The  relative  computational  size  of  the  problem  can  be  represented  by  the  total  number  of  grid  points  in 
a  model  volume  and  the  number  of  time  steps  in  the  model  run. 

Two  simulation  results  are  presented.  In  both  cases,  a  point  explosive  source  is  located  near  the  free 
surface.  The  first  example  is  for  a  plane  layered  geology.  The  second  example  uses  a  geology  that  contains  a  small 
topographic  feature,  penetrating  a  shallow  low-velocity  soil.  This  geology  is  common  in  a  number  of  problematic 
areas  of  the  world.  Both  results  show  complex  wave  phenomenon  that  can  be  exploited  for  development  of  system 
deployment  doctrines  or  in  system  adaptation.  Following  these  model  results,  we  discuss  preliminary  validation 
studies. 

3.1  Plane  Layered  Geology 

Figure  4  shows  two  y-z  plane  snapshots  for  the  vertical  component  of  particle  velocity.  The  model  used  an 
84  x  84  x  84  grid  volume  with  grid  separations  of  200  m.  The  center  frequency  of  the  excitation  source  was  roughly 
10  Hz.  In  Figure  4A  the  geologic  boundaries  are  indicated  by  horizontal  dotted  lines.  Table  1  gives  the  associated 
geologic  parameters.  This  model  result  is  representative  of  a  large  scale  crustal  velocity  structure.  Figure  4A  gives 
an  image  frame  at  0.008  s  into  the  simulation.  It  shows  that  the  explosive  pressure  point  source  has  been  initiated 
approximately  0.75  wavelengths  below  free  surface.  The  apparent  dipole  pattern  is  a  consequence  of  only  showing 
the  vertical  particle  motion.  Figure  4B  shows  a  snapshot  at  1.064  s.  A  variety  of  labeled  wavefronts  are  present.  The 
labeled  1  event  shows  an  upwardly  propagating  primary  reflection  from  L0  -  LI  layer  boundary,  2  shows  the 
primary  transmitted  wavefront  in  layer  L2.  Note  the  wavelength  elongation  in  this  higher  velocity  layer.  The  labeled 
3  event  shows  P-head waves  at  earth-air  interface  and  L0-L1  interface.  Labeled  wavefront  4)  shows  a  downwardly 
propagating  P-S  reflection  from  earth  air  interface.  And  barely  discemable  at  label  5  is  a  slight  amplitude  reduction 
on  the  primary  attributable  to  the  onset  of  the  boundary  absorption  condition. 


Layer  label 

Density  (kg/m3) 

V„(m/s) 

vs  (in/s) 

Grid  Dimension 

L0 

1500 

6000 

3400 

84  x  84  x  84 

LI 

2000 

7000 

4500 

L2 

2000 

8000 

5500 

Table  1.  Geologic  model  parameters  used  to  produce  the  results  shown  in  Figure  4. 
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Figure  4.  Vertical  slices  from  a  3D  model  result.  The  displays  show  the  vertical  component  of 
motion  resulting  from  an  explosive  point  source.  A)  Image  snapshot  0.008  s  after  the  explosive 
initiation.  Layer  interfaces  are  indicated  by  the  horizontal  dotted  lines  and  are  labeled  L0,  LI,  and 
Snapshot  image  at  1.064  s.  The  numbed  wavefronts  are  discussed  in  the  text. 
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3.2  Topographic  Intrusion 

A  large  percentage  of  the  world’s  most  troubled  regions  are  characterized  by  mountainous  terrain. 
Familiar  examples  include  North  and  South  Korea,  Taiwan,  the  Balkans,  Afghanistan,  and  northern  Iraq.  The  list  is 
extensive.  Independent  of  all  other  factors,  rugged  terrains  will  have  a  pronounced  effect  on  acoustic  and  seismic 
signal  properties.  One  of  the  unique  features  of  our  simulation  approach  is  the  ability  to  simulate  the  effects  of 
topography  on  seismic  waves.  In  Figure  5  we  illustrate  an  idealization  of  a  common  geologic  situation.  The  figure 
shows  a  rock  outcrop  penetrating  through  a  15-m-thick  compact  soil.  The  outcrop  extends  8  m  above  the 
surrounding  soil  horizon.  The  constituent  geologic  parameters  are  given  in  Table  2.  They  are  typical  of  values  seen 
at  shallow  depths.  The  explosive  source  is  placed  roughly  3  m  below  the  free  surface.  The  source  is  broadband  with 
a  center  frequency  of  roughly  40  Hz.  The  model  extent  is  165  x  156  x  62  m  and  uses  a  132  x  120  x  48  node  grid. 
There  were  2400  time  steps.  With  an  eight  processor  Sun  work  station  cluster,  the  runtime  was  roughly  7.5  hours. 


Layer  label 

Density  (kg/m3) 

V„(m/s) 

Vs  (m/s) 

Grid  Dimension 

L0 

1800 

800 

461 

132x120x48 

LI 

2400 

1800 

1040 

Table  2.  Geologic  an  model  parameters  used  to  produce  the  results  shown  in  Figure  6. 


An  animation  of  the  vertical  component  of  particle  motion  over  the  free  surface  and  in  a  cross  section  was 
generated  in  this  simulation.  Figure  6A  shows  a  snapshot  of  the  free-surface  at  0.097  s.  In  this  figure,  we  see  large- 
amplitude  surface  waves  radiating  from  the  source  point.  In  the  vicinity  of  the  topographic  feature,  we  see  that  faster 
traveling,  lower  amplitude  P-waves  have  refracted  through  the  hill.  These  waves  have  followed  a  path  of  shortest 
time  through  the  underlying  high  velocity  rock  layer.  Figure  6B  gives  a  snapshot  of  the  free-surface  at  roughly  0.2  s. 
In  this  time  frame,  we  clearly  see  that  the  hill  has  scattered  incident  waves  with  a  roughly  circular  radiation  pattern. 


On  the  lateral  (y)  sides  of  the  hill,  these  scattered  waves  are  interfering  with  surface  waves  that  are  diffracting 
around  the  hill.  This  complex  scattering  phenomenon  will  clearly  affect  sensor  performance. 

As  alluded  to  in  previous  sections,  and  as  seen  in  the  above  discussion,  topography  has  a  pronounced 
effect  on  seismic  signal  properties  and  will  therefore  affect  the  performance  of  systems  that  use  seismic  sensors.  To 
show  the  utility  of  simulations  in  predicting  system  performance,  we  use  the  topographic  simulation  results  to 
quantify  several  key  signal  characteristics  that  are  indicative  of  tracking  performance.  In  Figure  7  A  we  have  plotted 
the  vertical  component  of  ground  motion  as  a  function  of  time  for  a  60-element  line  array,  on  the  free-surface,  that 
spans  the  v  dimension  of  the  model.  Time  series  such  as  this  can  be  convolved  with  a  WAM  or  Raptor  sensor 
response  function.  The  resulting  data  can  then  be  used  in  WAM  or  Raptor  system  simulations  in  the  same  way  that 
field  data  are  used.  The  most  striking  characteristic  of  Figure  7 A  is  the  large  waveform  variation  across  the  sensor 
array.  The  waveforms  seen  in  the  vicinity  of  the  source  can  be  thought  of  as  the  “source  signature.”  In  the  vicinity 
of  the  hill  the  source  signature  has  been  radically  altered.  In  practical  terms  this  implies  that  target  features  needed 
for  classification  will  be  obscured.  In  Figure  7B  we  show  the  relative  peak-to-peak  amplitude  variation  across  the 
array.  Just  in  front  of  the  hill,  at  approximately  70  m,  there  is  an  abrupt  15  dB  amplitude  reduction  over  a  span  of 
less  than  7  m.  If  a  system  is  placed  in  this  region  and  it  uses  signal  level  to  estimate  target  range,  then  the  target  will 
appear  to  be  over  4  times  further  away  then  the  actual  range.  In  Figure  7C  we  plot  an  average  of  the  mean  signal 
coherence  (r)  at  the  source  center  frequency  (40  Hz).  In  this  coherence  plot,  a  moving  average  is  constructed  at  5-m 
intervals  using  a  20-m-long  subarray.  High  signal  coherence  is  needed  for  estimating  target  bearings.  In  the  plot  we 
see  that  the  coherence  function  is  generally  very  high.  However,  there  is  a  dramatic  drop  in  coherence  in  front  of  the 
hill.  In  this  region,  a  system’s  bearing  estimation  accuracy  will  be  severely  reduced. 

In  this  scenario,  all  three  major  ground  sensor  system  functions  will  be  adversely  affected  in  a  small  region 
in  front  of  the  hill.  An  analysis  approach  such  as  this  can  be  used  to  inform  the  user  of  where  optimal  system 
performance  can  be  achieved  for  any  given  geological  situation.  Furthermore,  system  developers  can  use 
information  such  as  this  to  select  robust  signal  characteristics  for  processing  and  to  specify  the  hardware  needed  to 
maintain  a  desired  performance  level. 


Figure  5.  A  common  geological  condition  in  mountainous  zones  is  a  hard-rock  intrusion  penetrating  a 
shallow  low- velocity  soil.  Model  results  are  obtained  for  the  idealized  geology  given  in  this  figure.  A)  and 
B)  show  oblique  and  cross-sectional  views  of  this  3D  scenario.  Table  2  gives  the  material  properties. 
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Figure  6.  A)  Vertical  component  of  ground  motion  image  at  0.097  s.  Concentric  surface  waves  are  radiating 
from  the  surface  point  directly  above  the  source.  Low  amplitude  P  waves  have  refracted  through  the  hill.  B) 
Image  frame  at  roughly  0.2  s.  Strong  surface  wave  diffraction  is  occurring  at  the  hill.  Low-amplitude  highly 
coherent  wavefronts  are  seen  moving  across  the  hill. 


Particle  velocity  w  (m/s);  h  =63.75  m;  <=0.207  s 


3.3  Preliminary  Numerical  Propagation  Validation 

The  FD-TD  propagation  model  must  accurately  represent  real  physical  processes.  An  extensive  field 
validation  is  planned  in  later  years  of  this  project.  Short  of  a  comparison  to  field  data,  preliminary  numerical 
validation  can  be  made  by  comparing  our  FD-TD  model  output  to  other  widely  accepted  numerical  results.  In  this 


discussion,  we  compare  seismograms  produced  by  our  program  to  seismograms  produced  by  the  wavenumber 
integration  program  OASP.  OASP  is  a  component  of  the  Oases  package  of  programs  (Version  2.1)  (Schmidt,  1997). 
It  uses  seismic  sources  in  horizontally  layered  media.  The  Oases  program  has  been  extensively  tested  and  results 
from  it  have  been  widely  reported  in  the  literature. 

Our  preliminary  OASES/FD-TD  comparison  uses  an  explosive  source  in  an  uniform  half-space.  The 
source  is  placed  at  a  depth  of  5  m  below  the  free  surface.  The  geologic  parameters  used  in  both  models  are:  P-wave 
velocity  of  700  m/s;  S-wave  velocity  of  400  m/s,  and  a  density  of  1800  kg/m3.  These  parameters  were  used  in  both 
the  FD-TD  and  OASP  models.  The  sensor  positions  are  shown  in  Figure  8 A.  An  overlay  of  the  vertical  particle 
velocity  records  for  the  surface  sensor  positions  are  shown  in  Figure  8B.  The  waveform  comparison  between  the 
two  methods  is  excellent. 

Another  check  on  the  FD-TD  model  results  can  be  made  by  comparing  the  P-to-S  arrival  time  curves  to 
ray  theory.  This  is  a  wave  that  leaves  the  source  as  a  P-wave,  hits  the  free  surface  at  an  oblique  angle,  and  is 
partially  converted  to  an  S-wave.  The  comparison  can  be  done  by  considering  a  snapshot  from  the  FD-TD  model 
result.  In  Figure  9A  we  show  a  snapshot  of  the  vertical  particle  velocity  at  time  0.1  s.  The  P-to-S  surface  reflection 
is  the  strong  linear  feature  that  runs  from  v  =  70  m  at  the  surface  to  v  =  30  m  at  a  depth  of  25  m.  Figure  9B  gives  a 
contour  plot  of  the  predicted  arrival  times  for  the  P-to-S  surface  reflection  based  on  ray  theory.  It  shows  the  same  P- 
S  linear  trend  as  that  observed  in  Figure  9A.  This  agreement  confirms  the  correctness  of  the  linear  wavefront 
character,  as  well  as  the  P-S  arrival  time. 


Figure  8.  A)  Snapshot  of  vertical  particle  velocity 
with  a  geophone  array  along  the  free  surface.  Note 
the  high  energy  Rayleigh  surface  wave  at  68  m.  B) 
Overlay  of  the  vertical  particle  velocity  records 
from  the  FD-TD  result  shown  in  A)  and  OASES 
waveforms.  The  first  pulse  on  each  record  is  the  P 
arrival,  and  the  second  is  the  Rayleigh  surface  wave. 
The  result  shows  excellent  agreement. 
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Figure  9.  A)  Image  snapshot  at  t=0.1s.  The  P-S  wave  is  the  high  amplitude  wavefront  shown  at  the  arrow. 
B)  P-S  arrival  time  contour  plot  drawn  from  ray  theory.  At  0.1  s,  ray  theory  predicts  that  the  P-S  wavefront 
will  located  at  the  position  shown  by  the  dotted  line.  This  agrees  well  with  the  location  of  the  P-S  wave 
seen  in  the  image  snapshot. 


4.0  Summary 

High  fidelity  simulated  data  can  reduce  system  development  costs  and  time.  In  this  paper  we  report  on  the 
mathematical  basis  for  a  new  FD-TD  method  for  simulating  propagating  seismic  signals.  The  technique  uses  a 
curved  coordinate  system  that  conforms  to  the  surface  topography  of  the  earth.  The  algorithm  has  been  parallelized 
following  a  domain  decomposition  strategy.  For  portability  it  is  written  in  Fortran77  and  uses  the  standard  MPI 
subroutine  library  to  distribute  sub-domains  across  a  processor  mesh.  The  code  has  been  successfully  run  on  clusters 
of  workstations  and  on  massively  parallel  Cray  and  IBM  platforms.  A  flexible  source  formulation  is  implemented 
that  allows  pressure  and  shear  force  excitations. 

Preliminary  results  are  given  for  a  plane  layered  model  and  for  a  rock-outcrop  penetrating  a  layer  of  soil. 
Both  models  show  complex  seismic  propagation  phenomena.  In  the  later  case,  topography  is  shown  to  have  a 
pronounced  effect  on  the  spatial  character  of  the  data.  In  particular,  waveforms  are  severely  altered,  amplitudes  are 
reduced  by  over  15  dB  over  distances  of  a  few  meters,  and  the  spatial  coherence  of  waveforms  is  reduced  from  0.95 
to  0.55  over  a  few  meters  distance.  All  three  of  these  characteristics  have  important  implications  for  ground  vehicle 
tracking  performance  in  battlefield  systems  that  use  seismic  data. 

A  preliminary  numerical  comparison  between  our  FD-TD  code  and  OASIS  wavenumber  integration 
models  shows  excellent  waveform  and  arrival  time  agreement.  Furthermore,  the  spatial  character  of  P-S  wavefronts 
and  the  P-S  wave  arrival  times  seen  in  our  FD-TD  model  are  in  excellent  agreement  with  ray  theory  predictions. 
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